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Efficient detection of a continuous-wave 
signal with a linear frequency drift. 
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ABSTRACT 

An efficient method i s presented for the detection of a continuous- wave 
(CW) signal with a frequency drift that is linear in time. Signals of this type 
occur if the transmitter and receiver are accelerating relative to one another, e.g. 
transmissions from the Voyager space craft. We assume that both the frequency 
and the drift are unknown. We also assume that the signal is weak compared to 
the Gaussian noise. The signal is partitioned into subsequences whose discrete 
Fourier transforms provide a sequence of instantaneous spectra at equal time 
intervals. These spectra are then accumulated* with a shift that is proportional to 
time. When the shift is equal to the frequency drift, the signal to noise ratio 
increases and detection occurs. In this paper we show how to compute these 
accumulations for many shifts in an efficient manner using a variant of the FFT, 
Computing time is proportional to LlogL where L is the length of the time 
series. Computational results are presented. 
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ABSTRACT 

An efficient method is presented for the detection of a continuous-wave 
(CW) signal with a frequency drift that is linear in time. Signals of this type 
occur in transmissions between any two locations that are accelerating relative to 
one another, e.g. transmissions from the Voyager space craft. We assume that 
both the frequency and the drift are unknown. We also assume that the signal is 
weak compared to the Gaussian noise. The signal is partitioned into subse- 
quences whose discrete Fourier transforms provide a sequence of instantaneous 
spectra at equal time intervals. These spectra are then accumulated with a shift 
that is proportional to time. When the shift is equal to the frequency drift, the 
signal to noise ratio increases and detection occurs. In this paper we show how to 
compute these accumulations for many shifts in an efficient manner using a vari- 
ant of the FFT. Computing time is proportional to LlogZ where L is the length 
of the time series. 
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1. Introduction 

For l=Q,...,L— 1 we are given a time series y t that is is generated from a weak 
continuous wave (CW) signal whose frequency is drifting linearly with time, i.e., 


S'/* 


2iri(u 0 + w,<i )/, 
At + 


( 1 - 1 ) 


The frequency at time t Q = 0 is co 0 and the frequency drift is Wj. Both are meas- 
ured in cycles per unit time. The term crg [ is Gaussian noise with mean zero and 
standard deviation a - . 

Signals of the form (1.1) result from CW transmissions between any locations 
that are accelerating with respect to one another, e.g. like those received from 
the Voyager space craft. They are also of interest for applications where the 
acceleration is considerably reduced but a high level of accuracy is required, such 
as satellite based aircraft navigational systems and similar systems that are 
currently under development for automobiles. Signals of this form are also 
thought to be the most probable form of initial communication from extraterres- 
trials [2]. 


THE PROBLEM: 

Given the time series y ; with a low signal to noise ratio then determine 
a> 0 and oij. 


Once a signal is detected, its verification and analysis is relatively simple using 
tuned versions of the algorithm that will be presented in the later sections. 
Therefore detection becomes the fundamental problem and the primary goal is to 
maximize its likehood. With a fixed computing resource we could persue the fol- 
lowing options: 

1. Develop a new method of detection. 

2. Speed the computation of an existing method and use higher resolution. 

3. Improve the accuracy of the existing method. 
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In this paper we will pursue options 2 and 3 in the context of an existing method 
which is described in [2]. They use the fast Fourier transform (FFT) to compute, 
but not to accumulate the spectra. Here we show that the FFT can be used in all 
of phases which results in a significant reduction in computing time. In particu- 
lar the FFT is used to compute certain spectral shifts with a high degree of accu- 
racy at the same time that it is used to accumulate the shifts. 

The approach in Cullers, et.al., [2] is to accumulate the spectra of subsequences 
to increase the signal to noise ratio. If L has factors L = MN we can define M 
subsequences y = y, where l = n + mN, for m = and 

n = 1. Each subsequence has N elements. Without loss of generality we 

can assume that the sequence y l is tabulated at integer times t f = l = n + mN 
and consequently (1.1) takes the form 

2'irt(o> 0 + mNiiiAmN 2‘ni(a> Q + 2mjVo), + n<i>,)n 

- Ac ' + (1-2) 

At this point, the sum of these subsequences would not increase the signal to 
noise ratio because y m n is not coherent as a function of m. Nevertheless, the fol- 
lowing steps can be taken: 

1. Compute the discrete Fourier transform (DFT) of the subsequences, i.e., 
compute M row transforms of length N 

2 TT 

1 N/2 — 1 — ikn 

= - 2 " ( 1 . 3 ) 

N n = - N/2 

If the frequency drift is negligible on each subsequence, i.e. if 
ncOj « o) fl + 2mNui v then Y k provides an approximation to the 
" instantaneous" spectrum at time mN. If we set ncOj = 0 in (1.2) and 
k = Af(o> 0 + 2mA r cu 1 ) in (1.3) then 


Y m,N{u 0 + 2mNu l ) ~ Ae 


2 tu ( oj 0 + mNdi^mN 


+ noise 


(1.4) 
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The contribution of the signal to Y m k is maximum at 
k = N(i i) Q + 2miVco 1 ) because the terms in (1.3) are coherent. 

2. Next compute the amplitude of the spectrum X . = Y . , and 

X m,N(*o + 2mNui x ) = 1 A 1 + n< ” 5C ( L 5) 


The purpose of this step is to provide an " instantaneous n spectra that is 
coherent for the next step. 

3. Finally the spectra X- . is summed with a "test" shift a ,i.e. 

J 


M— 1 

W- 2^ + „- 

m = 0 


( 1 . 6 ) 


From (1.5) and (1.6) the signal to noise ratio is maximum for 
k + am = AT(co 0 + 2mATa) 1 ) or for k = Niii Q and a = 2N 2 (x)y Therefore the 
problem becomes one of computing the k and a that correspond to the max- 
imum of £ A (a) since a) Q is then given by oj 0 = k/N and (Oj is given by 
Wj = a/(2AT ). To implement (1.6) it is necessary to interpolate the spectra to 
non- integer values k + am of the subscript. An accurate method for this interpo- 
lation is given in the next section. 

The maximum of (1.6) is located by tabulating ^(a) at a set of points 
oty = a + j 8 for Accuracy can be improved by a second tabulation 

in the neighborhood of the maximum of the first tabulation. This would establish 
verification and improve the accuracy of the analysis. If we set S -. = S.( a ) 
then in the sections that follow we will be concerned with the efficient tabulation 
of a discrete version of (1.6) namely (2.1). 


2. The Algorithm. 
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Given a, 8 and the MxN matrix X m n , then we wish to compute another MxN 
matrix 

M— 1 

Sj,k ~ 2 ■^m,Jfc + m(a+/ 8 )‘ ( 2 - 1 ) 

m = 0 


The jth row of S- . contains the sum of the rows of X . where the mth row is 

j m,K 

shifted by m(a+j'8). The smallest shift between any two adjacent rows is 8. 
Fractional shifts are permitted and computed with a high degree of accuracy 
using trigonometric interpolation and the FFT. 


Let x be the rowwise backward Fourier transform of X , . Then 

n tTljK 




AT/ 2-1 
y x t 

xLi m,n 

n = — N/2 


— ikn~ 


2tt 

N 


( 2 . 2 ) 


and 


L m,Jt + m(a+/8) 


N 


N/ 2-1 

2 

n = -N / 2 


— t [A + m (a+ j 8 )] rr 


X m t 
m } n 


2tt 

N 


(2.3) 


From (2.1) 




N 


N/ 2-1 

s 

n = - N/2 


2tt 

M—l — imno. 

N 




2 ( c 

m = 0 


„) C 

in, n * 


2tt 

N 


— ikn~ 


2it 

N 


(2.4) 


or 



N 


N/ 2-1 

2 

n = — N/2 


— ikn~ 


Jt n 


2tt 

N 


(2.5) 


where 


M-l 


— 18 /mrr 




2tt 


m =0 


( 2 . 6 ) 


and 
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— imna" 


Z m,n = C 


2tt 

N 


m,n 


(2.7) 


Given a, 8, and X m k , the algorithm consists of the following steps: 

A. Compute x m n from the inverse of (2.2) and z m n from (2.7). 

B. Next compute «. using (2.6). With fixed n, each column can be 
transformed independently using a method that will be presented below. 

C. Compute S- k from (2.5) using the forward FFT applied to each row. 

D. If Sj max ^ max = max S- k then o) Q = kmax/N and co 1 = (a + jmaxh)/(2N 2 ) 

3 


The FFT can be used for steps A and C. However it cannot be used for step B 
since the complex exponential contains 8 which is not necessarily an integer. 
Nevertheless we will show that (2.6) can be computed in about the time required 
by a fast AZ-point convolution. We will use a technique introduced by Bluestein 
[1] that enabled him to develop an FFT for arbitrary M including large prime 
numbers. The details of his algorithm are presented in [4] as well as its imple- 
mentation on multiprocessors. 

The same transform (2.6) is applied to each column of the array z m n and hence 
the exposition can be simplified by removing the subscript n. Therefore we wish 
to compute 

2it 

M— 1 — t (3 jm 

N (2-8) 

m = 0 

where fj = 8n. We begin with the following identity 

-2 jm = (j-mf - j 2 - to 2 . (2.9) 


Substituting into (2.7) we obtain 
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( 2 . 10 ) 

( 2 . 11 ) 

( 2 . 12 ) 

(2.13) 


Equation (2.11) is preferable to (2.8) for computing Sy because (2.11) can be 
evaluated efficiently using the FFT. However an indirect approach is required 
because (2.11) is not a circular convolution, i.e., the subscript j—m is not inter- 
preted modulo M. If j—m is negative then it is replaced by m—j rather than 
j—m + M. Nevertheless, one can construct an equivalent circular convolution. 

Select P^2M—2 as a highly factorizable integer, e.g., a power of two. Next, 
extend b m and c m to sequences with P elements as follows: 


3 

II 

o 

M^m <P — M 

(2.15a) 

b = t p m 

T7i r — m 

P — m <P. 

(2.15b) 

c* 

2 

II 

O 

M< m <P 

(2.15c) 

Define 



A 

a i = 

P- 1 

= ft." 1 y 6. c 

j jLj j — m m 

m = Q 

(2.16) 

where the subscripts are now interpreted modulo P. It can be determined by 
inspection that a ■ = a . for j = 1. Since (2.16) is a P point circular 

convolution, it can be evaluated using the familiar FFT procedure [4]. The algo- 
rithm for fast evaluation of (2.11) can now be summarized. 
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B.O Compute b m from (2.12) and extend to length P using (2.15a) and (2.15b). 
Perform a P point forward FFT and call the result B k . 

B.l Compute c m from (2.13) and extend to length P using (2.15c). Perform a P 
point forward FFT and call the result C k . 

B.2 Compute D k — B k C k . Perform a backward P-point FFT and call the result 

<v 

B.3 Multiply the first M elements of d. by b- 1 to obtain a- = a-. The remain- 

* J J J 

ing P — M elements of d k may be discarded. 


Note that step B.O is an initialization step that does not have to be repeated for 
the analysis of any subsequent time series. 


3. Computational Tests 

The algorithm described above has been implemented for test purposes. The 
resulting program is written entirely in Fortran without assembly code or 
assembly-coded library routines. Some modest effort was made to insure that the 
code was suitable for reasonably high performance execution on vector comput- 
ers, but in general it is an entirely straightforward implementation of the algo- 
rithm presented in section 2. 

Two details of this implementation are worth mentioning because they would be 
worth inclusion in any serious production system. First of all, every FFT opera- 
tion, including those that are a part of the fractional FFT computation, was per- 
formed using “simultaneous” FFT routines. In other words, the FFTs were per- 
formed vector-wise instead of individually. Such an implementation is of course 
highly suited for a vector computer, but it also is quite appropriate for imple- 
mentation on a parallel computer system. Secondly, it is important to note that 
a substantial amount of computation time can saved by employing real-to- 
complex FFTs in step A above and by employing complex-to-real FFTs in step C 
above. A by-product of this modification is that the FFTs required for the frac- 
tional transforms in step B need only be performed on iV/2+1 columns. 
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The detection of the true signal base frequency and drift rate is enhanced in this 
implementation by employing a weighted score. This is because even when a 
pure signal with no noise is input to the processing algorithm, the peak in the 
result array is “smeared” over a number of nearby elements. This pattern of 
smearing from a pure signal can in fact be used to design an improved detection 
filter. The authors found that the following weighted score formula proved both 
inexpensive and effective: 


T j,k 1 .0525 ; _ 2k + 1 + S ; - _ j ^ k + Sj _ lfk + 1 + 1.2WS jtk + S. +hk _ 1 +S. +hk (3.1) 


+ 1.0525 


j + 2,k-l 


where 5 - ^ is the final result of step C above. The statistical significances of 
these weighted scores are determined by explicitly computing the mean and stan- 
dard deviation of all T ■ , . 

Table 1 contains a detailed accounting of the computational cost of the complete 
detection algorithm. “C-C”, “R-C” and “C-R” denote the three types of FFTs: 
complex-to-complex, real-to-complex, and complex-to-real. The column headed 
“Ref.” lists references to specific equations and algorithm steps. The colu mn 
headed “Operation Counts” contains the number of real floating point arithmetic 
operations for each step, with adds, subtracts and multiplies each counting as one 
operation. Square roots, which are required to compute amplitudes, are counted 
as 15 floating point operations (this approximates the cost of a typical 64-bit 
square root function). 
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Table 1 

Floating point operation counts 

Computational Step 

Ref. 

Operation Count 

Initial row-wise C-C FFTs 

(1.3) 

5MN\og 2 N 

Amplitudes 

(1.5) 

1SMN 

Inverse row- wise R-C FFTs 

A 

{5MN\og 2 N)/2+MN 

Column-wise multiplications 

B.l 

3MN 

Column- wise forward C-C FFTs 

B.l 

5MN\og 2 M+5MN 

Pointwise multiplications 

B.2 

6MN 

Column- wise inverse C-C FFTs 

B.2 

5MN\og 2 M+5MN 

Column-wise multiplications 

B.3 

3MN 

Forward row- wise C-R FFTs 

C 

(5M/Vlog 2 7V)/2 + MN 

Weighted scores and statistics 
Total 

(3.1) 

13MN 

MN(10\og 2 MN+55) 


This algorithm was exercised by performing a series of tests with pseudoran- 
domly selected initial frequencies and drift rates. In these tests, M and N were 
each set to 1,024. The initial value of the drift rate a was set to — 1/2, and the 
drift rate increment, 8, was set to -1/M. Thus the sums S k ( a) in (1.6) were 
computed for all drift rates a from —1/2 to 1/2. For each test, a linearly drifting 
signal was added to complex Gaussian pseudorandom data of the form 
9 1 — ( e / + / where e t and f l are real Gaussian data with zero mean and 

unit variance. The ratio of the amplitudes of the signal and the noise was 56.6, 
so that the power of the signal was lower than that of the noise by a ratio of 
3,200. This extremely low signal to noise level is approximately the level that the 
Cyclops SETI system hopes to detect [2], 
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Table 2 

Algorithm Test Results 

Trial 

Generated 

Detected 



No. 

Base Freq. 

Drift Rate 

Base Freq. 

Drift Rate 

if -score 

Time 

1 

119.842 

-0.26187 

119.00 

-0.2607 

6.506 

2.3448 

2 

967.298 

0.49238 

968.00 

0.4922 

6.515 

2.4654 

3 

706.734 

0.04648 

706.00 

0.0479 

6.570 

2.4522 

4 

828.500 

0.05203 

829.00 

0.0518 

7.146 

2.4778 

5 

97.555 

0.35731 

98.00 

0.3564 

6.991 

2.5006 

6 

626.842 

-0.25255 

627.00 

-0.2529 

6.303 

2.5685 

7 

989.664 

-0.14073 

989.00 

-0.1396 

7.688 

2.5031 

8 

203.058 

-0.42424 

203.00 

-0.4248 

8.274 

2.3588 

9 

1022.170 

-0.10787 

1022.00 

-0.1074 

7.253 

2.4729 

10 

921.632 

0.43326 

922.00 

0.4336 

8.863 

2.4977 

Ave. 





7.211 

2.4642 


These tests were performed on the Cray-2 operated by the NAS Systems Division 
at NASA Ames. The results of these tests are displayed in Table 2. The column 
headed “Z-score” gives the Z-score of the detection, i.e. the number of standard 
deviations above the mean. The column headed “Time” contains single proces- 
sor CPU times in seconds. These results indicate that the detection algorithm is 
effective in determining the unknown base frequency and drift rates. In every 
trial, the element of the final array with the highest weighted score was within 
one bin both in drift rate and frequency from the true value. These Z-scores are 
well above the level of random scores since one would not expect a random Z- 
score to exceed 5 in a 1,024X1,024 array. 

The average processing CPU time, 2.464 seconds, corresponds to a performance 
rate of approximately 109 MFLOPS, based on the total number of floating point 
operations given in Table 1. This performance rate could be increased by 
employing assembly-coded simultaneous FFT routines and by employing all four 
of the Cray-2 processors. In any event, these timings indicate that such 
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computations are feasible given suitably fast computer hardware. 

As was mentioned above, the central problem is making an accurate initial detec- 
tion. The algorithm presented here accomplishes this basic objective. Once a 
putative detection has been made, its confidence level can be increased in a 
number of ways, such as by repeating our algorithm with higher resolution near 
the putative detection, or by employing a matched filter keyed to the putative 
frequency and drift rate. It is also expected that some modest enhancements and 
improvements in the algorithm could also improve these detection levels. 
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